data = readtable('../../data/main_VAR/UK_source_pre1946.xlsx', 'Sheet', 'Full_Data_for_VAR');

date = data.Year;
date0 = date;

startdatenum = 1729;
enddatenum = 1914;
middatenum = 1794;

startdate = find(date == startdatenum);
middate = find(date == middatenum);
enddate = find(date == enddatenum);
T_post = length(date(startdate:enddate));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Load data
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

option.figure = 0;

Tstart = startdate; 
Tend = enddate; 
T = Tend - Tstart + 1;

date = date(startdate:enddate);

data = data(find(date0 >= startdatenum - 1 & date0 <= enddatenum), :);

taxrevgdp = data.revtogdp(2:end); % government tax revenue to GDP ratio
taxrevgdp0 = data.revtogdp(1); % government tax revenue to GDP ratio
spendgdp = data.spendingtogdp(2:end); % government spending before interest exdp. to GDP ratio
spendgdp0 = data.spendingtogdp(1); % government spending before interest exp. to GDP ratio
surplusgdp = taxrevgdp - spendgdp; % primary surplus to gdp
rpcgdpgr = data.realGdpGrowth(2:end); % real gdp growth (in logs)
inflation = data.inflation(2:end); % growth in GDP price deflator (in logs)

ynom1 = data.short_r(2:end) / 100; % nominal yield on a 1-year Treasury bill, expressed per annum
ynom10 = data.rate_10year(2:end) / 100; % nominal yield on a 10-year Treasury bond, expressed per annum
dy = data.dy(2:end);
pdm = log(1 ./ dy); % log price/dividend ratio

%% convenience yield

cy_short1 = xlsread('../../data/supplement/UK_cony.xlsx', 'UK basis (ST)', 'ck5:ck45'); % 1873-1913
cy_short2 = xlsread('../../data/supplement/UK_cony.xlsx', 'UK basis (ST)', 'ck57:ck63'); % 1925-1931

cy_short = [mean(cy_short1) * ones(1873 - startdatenum + 1, 1); cy_short1] / 100;
cy_short = [cy_short; mean(cy_short)];

cy_long1 = xlsread('../../data/supplement/UK_cony.xlsx', 'UK basis (LT)', 'ck5:ck45'); % 1873-1913
cy_long2 = xlsread('../../data/supplement/UK_cony.xlsx', 'UK basis (LT)', 'ck57:ck63'); % 1925-1931

cy_long = [mean(cy_long1) * ones(1873 - startdatenum + 1, 1); cy_long1] / 100;
cy_long = [cy_long; mean(cy_long)];

% cy is zero before 1794;
cy_long = [zeros(middate - startdate + 2, 1); cy_long(3 + middate - startdate:end)];
cy = cy_long;

ynom1 = ynom1 + [zeros(middate - startdate + 1, 1); cy_short(3 + middate - startdate:end)];
ynom10 = ynom10 +cy_long(2:end);

ynom1 = log(ynom1 +1);
ynom10 = log(ynom10 +1);

divgrm = diff(log(data.dy) + data.stockindex);

gdebt = data.debttogdplevel;

% calibrate cy/gdp ratio
cy_data = mean((cy_long1) / 100 .* gdebt(1873 - startdatenum + 2:1913 - startdatenum + 2));

datew0 = (1728:1914)';
seig_rev = zeros(size(datew0));
seig_rev(1:find(datew0 == 1794)) = 0;
seig_rev(find(datew0 == 1795):find(datew0 == 1872)) = cy_data;
seig_rev(find(datew0 == 1873):find(datew0 == 1913)) = gdebt(find(datew0 == 1873):find(datew0 == 1913)) .* (cy_long1 / 100);
seig_rev(find(datew0 == 1914)) = cy_data;

% adjust for convenience yield
taxrevgdp = taxrevgdp + seig_rev(2:end);
taxrevgdp0 = taxrevgdp0 + 0; % zero convenience yield in 1728

deltalogg = [log(spendgdp(1)) - log(spendgdp0); log(spendgdp(2:end)) - log(spendgdp(1:end - 1))];
deltalogtau = [log(taxrevgdp(1)) - log(taxrevgdp0); log(taxrevgdp(2:end)) - log(taxrevgdp(1:end - 1))];

yieldspr = ynom10 - ynom1; % spread between 10-yr and 1-yr yields
y0nom_1 = mean(ynom1);
y0nom_10 = mean(ynom10);
yspr0 = mean(yieldspr);
pi0 = mean(inflation);
x0 = mean(rpcgdpgr);
tau0 = mean(deltalogtau);
g0 = mean(deltalogg);
surplus0 = mean(surplusgdp);

% log div/gdp ratio growth
divgrm = divgrm - inflation;

A0m = nanmean(pdm);
% mu_m = nanmean(divgrm);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Detrend tax and spending
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

tts(1) = taxrevgdp(1); 
gts(1) = spendgdp(1); 

for t = 2:T
    gts(t) = gts(t - 1) * exp(deltalogg(t) - g0);
    tts(t) = tts(t - 1) * exp(deltalogtau(t) - tau0);
end


run ../../tools/run_var_estimate;


save('MAT/res_annual_VAR_nodebt_cy.mat','-regexp', '^(?!(option|Globaloption)$).');
